Systematic study of deformed nuclei at the drip lines and beyond 
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An improved prescription for choosing a transformed harmonic oscillator (THO) basis for use in 
configuration-space Hartree-Fock-Bogoliubov (HFB) calculations is presented. The new HFB+THO 
framework that follows accurately reproduces the results of coordinate-space HFB calculations for 
spherical nuclei, including those that are weakly bound. Furthermore, it is fully automated, facili- 
tating its use in systematic investigations of large sets of nuclei throughout the periodic table. As 
a first application, we have carried out calculations using the Skyrme Force SLy4 and volume pair- 
ing, with exact particle number projection following application of the Lipkin-Nogami prescription. 
Calculations were performed for all even-even nuclei from the proton drip line to the neutron drip 
line having proton numbers Z = 2,4,..., 108 and neutron numbers TV = 2,4,..., 188. We focus 
on nuclei near the neutron drip line and find that there exist numerous particle-bound even-even 
nuclei (i.e., nuclei with negative Fermi energies) that have at the same time negative two-neutron 
separation energies. This phenomenon, which was earlier noted for light nuclei, is attributed to 
bound shape isomers beyond the drip line. 
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I. INTRODUCTION 

The development of experimental facilities that accel- 
erate radioactive ion beams Jj,,^] has opened up a window 
to many nuclei that were heretofore inaccessible. With 
these new facilities and the new detector technology that 
is accompanying them, it is becoming possible to study 
the properties of nuclei very far from the valley of beta 
stability, all the way out to the particle "drip lines" and 
perhaps even beyond. 

Much work is now in progress to develop appropri- 
ate theoretical tools for describing nuclei in these exotic 
regimes 0]. A proper theoretical description of such 
weakly-bound systems requires a careful treatment of 
the asymptotic part of the nucleonic density. An appro- 
priate framework for these calculations is Hartree-Fock- 
Bogoliubov theory, solved in coordinate representation 
P, 1^, 1^ . This method has been used extensively in the 
treatment of spherical systems but is much more difh- 
cult to implement for systems with deformed equilibrium 
shapes (7|,^, 

In the absence of reliable coordinate-space solutions 
to the deformed HFB equations, it is useful to consider 
instead the configuration-space approach, whereby the 
HFB solution is expanded in a single-particle basis. One 
approach has been to use a truncated basis composed 
partly of discrete localized states and par tly of discretized 
continuum and oscillating states 0,11,^3 • Because of the 
technical difficulties in implementing this method, it has 
typically been restricted to include states in the contin- 
uum up to at most several MeV. As a consequence, such 
an approach should not be able to describe adequately 



the spatial properties of nuclear densities at large dis- 
tances. 

An alternative possibility is to expand in a basis of 
spatially localized states. Expansion in a harmonic os- 
cillator (HO) basis is particularly attractive because of 
the simple properties of oscillator states. There have 
been many configuration-space HFB+HO calculations 
reported, either employing Sky rme forces or the Gogny 
effective interaction [ll|, HJI ■, or using a relativistic 
Lagrangian 0, ^ . This methodology has proven par- 
ticularly useful when treating nuclei in or near the valley 
of stability. For nuclei at the drip lines, however, the 
HFB+HO expansion converges slowly as a function of 
the number of oscillator shells [6| , producing wave func- 
tions that decrease too steeply at large distances. The 
resulting densities, especially in the pairing channel, are 
artificially reduced in the outer region and do not reflect 
correctly the pairing correlations of these weakly-bound 
nuclei. 

A related approach that has recently been proposed 
is to instead expand the quasiparticle HFB wave func- 
tions in a complete set of transformed harmonic oscil- 
lator (THO) basis states [Ulil 111, obtained by ap- 
plying a local-scaling coordinate transformation (LST) 
111 0, i2 to the standard HO basis. The THO basis 
preserves many useful properties of the HO wave func- 
tions, including its simplicity in numerical algorithms, 
while at the same time permitting us to incorporate the 
appropriate asymptotic behavior of nuclear densities. 

Applications of this new HFB+THO methodol ogy 
have been reported both in the non-relativistic [T^ . Il9| 
and relativistic domains In all of these calculations. 
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specific global parameterizations were employed for the 
scalar LST function that defines the THO basis. There 
are several limitations in such an approach, however. On 
the one hand, any global parameterization of the LST 
function will of necessity modify properties throughout 
the entire nuclear volume, in order to improve the asymp- 
totic density at large distances. This is not desirable, 
however, since the HFB+HO results are usually reliable 
in the nuclear interior, even for weakly-bound systems. 
In addition, because of the need to introduce matching 
conditions between the interior and exterior regions, a 
global LST function will invariably have a very compli- 
cated behavior, especially around the classical turning 
point, making it difficult to simply parameterize. Per- 
haps most importantly, the minimization procedure that 
is needed in such an approach to optimally define the ba- 
sis parameters is computationally very time consuming, 
especially when a large number of shells is included, mak- 
ing it very difficult to apply the method systematically 
to nuclei across the periodic table. 

In the present work, we propose a new prescription for 
choosing the THO basis. For a given nucleus, our new 
prescription requires as input the results from a relatively 
simple HFB-I-IIO calculation, with no variational opti- 
mization. The resulting THO basis leads to HFB+THO 
results that almost exactly reproduce the coordinate- 
space HFB results for spherical Q and axially deformed 
[Toll nuclei and are of comparable quality to those of the 
former, more complex, HFB-f THO methodology . 

Because the new prescription requires no variational 
optimization of the LST function, it can be readily ap- 
plied in systematic studies of nuclear properties. As the 
first such application, we carry out a detailed study of 
nuclei between the two-particle drip lines throughout the 
periodic table, using the Skyrme force SLy4 [2^ and 
volume pairing 0. In order to restore good particle 
number, we apply the Lipkin-Nogami (LN) prescription 
mm 111 113^ [III followed by exact particle- number 
projection (PNP) |30| . 

The structure of the paper is the following. In Sec. 
IhI we briefly review the HFB and LN methods, noting 
several features particular to its coordinate and config- 
urational representation. In Sec. IIIII we introduce the 
THO basis and then formulate our new prescription for 
the LST function. The results of systematic calculations 
of even-even nuclei are reported in Sec. IIVI with special 
emphasis on those nuclei that are at the neutron drip 
line and just beyond. Conclusions and thoughts for the 
future are presented in Sec. 



II. OVERVIEW OF 
HARTREE-FOCK-BOGOLIUBOV THEORY AND 
THE LIPKIN-NOGAMI METHOD 

In this section, we review the basic ingredients of 
Hartree-Fock-Bogoliubov theory and the Lipkin-Nogami 
method followed by particle-number projection. Since 



these are by now standard tools in nuclear structure, 
we k eep the presentation brief and refer the reader to 
Ref. [301 for further details. 

HFB is a variational theory that treats in a unified 
fashion mean-field and pairing correlations. The HFB 
equations can be written in matrix form as 

where _B„ are the quasiparticle energies, A is the chemical 
potential, h = t + T and A are the Hartree-Fock (HF) 
hamiltonian and the pairing potential, respectively, and 
Un and Vn are the upper and lower components of the 
quasiparticle wave functions. These equations are solved 
subject to constraints on the average numbers of neu- 
trons and protons in the system, which determine the 
two corresponding chemical potentials, A„ and Ap. 

In coordinate representation, the HFB approach con- 
sists of solving (|2.1|) as a set of integro-differential 
equations with respect to the amplitudes [/(£"„, r) and 
V(En,v), both of which are functions of the position co- 
ordinate r. The resulting density matrix and pairing ten- 
sor then read 

p(r,r') = ^*(^n,r)y(i?„,r') , (2.2a) 

<£„<£,„ ax 

K(r,r') = J2 V*iEn,r)UiEn,v') . i2.2h) 

Typically, the HFB continuum is discretized in this ap- 
proach by putting the system in a large box with appro- 
priate boundary conditions 0. 

In the configurational approach, the HFB equations 
are solved by matrix diagonalization within a chosen 
single-particle basis {ipa} with appropriate symmetry 
properties. In this sense, the amplitudes Un and Vn enter- 
ing Eq. (|2.1|) may be thought of as expansion coefficients 
for the quasiparticle states in the assumed basis. The nu- 
clear characteristics of interest are determined from the 
density matrix and pairing tensor, 

p(r,r') = ^p„^Va(r)^^(r') , (2.3a) 

a/3 

K(r, r') = YKa0ipa{r)^i3{r') , (2.3b) 

a(3 

which are expressed in terms of the basis states tpa and 
the associated basis matrix elements as 

Po.p= E V:JEn)V0n{En) , (2.4a) 
o<B„<£;„ax 

= E V*n{E^)Ui3n{E„) . (2.4b) 

o<£;„<£;,„ax 

In configuration-space calculations, all quasiparticle 
states have discrete energies En- 
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The results from configuration-space HFB calculations 
should be identical to those from the coordinate-space 
approach when all the states ipa from a complete single- 
particle basis are taken into account. Of course, this is 
never possible. In the presence of truncation, it is essen- 
tial that the basis produce rapid convergence, so that reli- 
able results can be obtained within computational limita- 
tions on the number of basis states that can be included. 

The LN method serves as an efficient method for 
restoring particle number before variation ,241] . With 
only a slight modification of the HFB procedure outlined 
above, it is possible to obtain a very good approxima- 
tion for the optimal HFB state, on which exact par ticle 
number projection then has to be performed |28l l3lj|. 

In more detail, the LN method is implemented by per- 
forming the HFB calculations with an additional term 
included in the HF hamiltonian, 



h' ^h- 2X2(1 -2p), 



(2.5) 



and by iteratively calculating the constant A2 (separately 
for neutrons and protons) so as to properly describe the 
curvature of the total energy as function of particle num- 
ber. For an arbitrary two-body interaction V, A2 can be 
calculated from the particle-number dispersion according 

to il. 



A. 



(0|T/|4)(4|jV^|0) 
(0|7V2|4)(4|iV2|0) 



(2.6) 



where |0) is the quasiparticle vacuum, N is the parti- 
cle number operator, and |4)(4| is the projection oper- 
ator onto the 4-quasiparticle space. On evaluating all 
required matrix elements, one obtains p7| 



A2 



4Trr'/9(l - p) + 4TrA'(l - p)k 



8 [Trp(l - p)Y 
where the potentials 



16Trp2(i_p)2 



r' 



A' 



E 



V^„f,'„'{p{l - p))u'u. 



(2.7) 

(2.8a) 
(2.8b) 



can be calculated in full analogy to F and A by replacing 
the p and k in terms of which they are defined by p(l — p) 
and respectively. In the case of the seniority pairing 
interaction with strength G, Eq. (|2.7(l simplifies to 



A. 



G Tr(l - p)k TrpK - 2 Tr(l 



4 [Trp(l-p)]^-2Trp2(l^p)2 



(2.9) 



An explicit calculation of A2 from Eq. (|2.7|l requires 
calculating new sets of fields Eq. (|2.8|l . which is rather 
cumbersome. However, we have found that Eq. (|2.7|) 
can be well approximated by the seniority-pairing expres- 
sion Eq. (|2.9|l with the effective strength 



G — Gc 



A2 



(2.10) 



determined from the pairing energy 
^^pair = -^TrAK 
and the average pairing gap 



A = 



TrAp 



(2.11) 



(2.12) 



The use of the LN method in HFB theory requires spe- 
cial consideration of the asymptotic properties of quasi- 
particle states 0) ) of essential importance for weakly- 
bound systems. Because of the modified HF hamilto- 
nian 12.5|l . new terms appear in the HFB-I-LN equation, 
which are non-local in coordinate representation and thus 
can modify the asymptotic conditions. Effectively, this 
means that the standard Fermi energy A has to be re- 
placed by 



or by 



A' = A + 2A2(l-2nmin) 



A" = A + 2A2 



(2.13) 



(2.14) 



pair 



where nmin is the norm of the lower HFB component 
^(-Einin, »") corresponding to the smallest quasiparticle 
energy £'i„i„. 

The first expression 1)2.13(1 assumes that the asymp- 
totic properties can be inferred from the HFB equation 
in the canonical basis, in which p is diagonal and has 
eigenvalues that can be estimated by norms of the sec- 
ond HFB components. The second expression (|2.14|) per- 
tains to the HFB equation in coordinate representation, 
in which the integral kernel p(r^ r') vanishes at large dis- 
tances. Neither of these expressions can be rigorously jus- 
tified, thereby demonstrating limitations of using the LN 
method to analyze spatial properties of wave functions. 
These ambiguities are enhanced by the fact that the LN 
method overestimates the curvature A2 near magic num- 
bers miiii. 

Note that in the exact projection before variation 
method, the Fermi energy is entirely irrelevant, and hence 
one should not attribute too much importance to the 
choice between A' and A". Nevertheless, since the PNP 
affects only occupation numbers, leaving the canonical 
wave functions unchanged, in what follows we use the 
modified Fermi energy A' in modelling the asymptotic 
behavior needed to implement the THO method. 

Finally, we should note that the HFB machinery 
detailed above can be readily implemented with a 
quadrupole constraint [33| , as is the case for some of the 
calculations we will be reporting. 



III. THE TRANSFORMED HARMONIC 
OSCILLATOR BASIS 

In the present study, we carry out HFB calculations 
in configuration space, expanding in a transformed har- 
monic oscillator basis. This basis was originally intro- 
duced in Refs. fvA ITM IT^ . and we refer the reader to 
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Ref. 19*1 for details concerning the use of the deformed 
THO basis and for a discussion of the cut-off procedure 
that is used to perform the summations in Eq. p. 4(1 . We 
also refer the reader to an interesting new application of 
the THO basis to one-dimensional problems of interest 
in molecular physics [s^l- 

As noted earlier, all previous calculations using the 
THO basis in HFB calculations employed a global pa- 
rameterization of the LST function that defined the ba- 
sis. In the following subsections, we develop a new and 
improved form for the transformation, which we then use 
in the HFB+THO applications to be reported in Sect. 



A. Comparison of coordinate-space HFB 
calculations and configuration-space HFB-|-HO 
calculations 

The main differences between the results of coordinate- 
space HFB calculations and those from configuration- 
space HFB-I-HO calculations can be seen in plots of the 
corresponding local density distributions. A typical ex- 
ample is shown in Fig. ^ where the densities and their 
logarithmic derivatives from coordinate-space HFB cal- 
culations (solid lines) are compared with those from a 
configurational HFB-I-HO calculation. Although the cal- 
culations were done for a specific spherical nucleus and 
Skyrme interaction, the features exhibited are generic. 
Note that the coordinate-space HFB calculations were 
carried out in a box of 30 fm, so that the logarithmic 
derivative of the density obtained in that calculation 
shows a sudden drop near the box edge. 

Invariably, the logarithmic derivative p' j p associated 
with the coordinate-space HFB solution shows a well- 
defined minimum near some point -Rmin in the asymptotic 
region, after which it smoothly approaches a constant 
value —A;, where 

k^2K^ 2\j2m(E^,^ - \')/h^ (3.1) 

is associated with the HFB asymptotic behavior for the 
lowest quasiparticle state that has the corresponding 
quasiparticle energy E'min (see Eq. H2.13|l and Ref. 0). 
This property is clearly seen in the upper panel of Fig. ^ 
One can also see that the HFB-I-HO densities and loga- 
rithmic derivatives are in almost perfect agreement with 
the coordinate-space results up to (or around) the dis- 
tance i?min- We conclude, therefore, that the HFB+HO 
densities are numerically reliable up to that point. 

Moreover, the HFB value of the density decay constant 
k=2K, when calculated from Eq. (|3.1|l . is also correctly 
reproduced by the HFB-I-HO results. It is not possible 
to distinguish between the values of k that emerge from 
the coordinate-space and harmonic-oscillator HFB calcu- 
lations, both values being shown by the same line in the 
upper panel of Fig. ^ 

Soon beyond the point i?min, the HFB-I-HO density 
begins to deviate dramatically from that obtained in the 
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FIG. 1: Logarithmic derivative of the density (upper panel), 
and the density in logarithmic scale (lower panel), as func- 
tions of the radial distance. The coordinate-space HFB re- 
sults (solid line) are compared with those for the HFB-I-HO 
method (denoted p) with A'sh^ 8, 12 and 20 HO shells, as 
well with the approximation (denoted p) given by Eq. 13.711 
(small circles). 



coordinate-space calculation. For relatively small num- 
bers of harmonic oscillator shells A'gh, the logarithmic 
derivative of the HFB-I-HO density goes asymptotically 
to zero following the gaussian behavior of the harmonic 
oscillator basis. The resulting HFB-I-HO density does 
not develop a minimum around the point -Rmin, as seen 
from the A'sh = 8 curve shown in the upper panel of 
Fig. When the number of harmonic oscillator shells 
A'sh increases, the HFB-I-HO solution tries to capture 
the correct density asymptotics. Due to the gaussian 
asymptotic of the basis, however, the logarithmic deriva- 
tive of the HFB-I-HO density only develops oscillations 
around the exact solution (see the A^sh = 12 and 20 
curves in the upper panel of Fig. As a result, the 
logarithmic derivative of the HFB-f HO density is very 
close to the coordinate-space result around the mid point 
Rm = (-Rmax — -Rmin)/2, whcrc i^max IS the positiou of the 
first maximum of the logarithmic derivative after i?min- 

In summary, the following HFB+HO quantities agree 
with the coordinate-space HFB results: (i) the value of 
the density decay constant fc; (ii) the local density up 



5 



to the point i?niin where the logarithmic derivative p' j p 
shows a clearly- defined minimum; (in) the actual value of 
this point i?min; i'>''v) the value of the logarithmic deriva- 
tive of the density at the point i?m defined above. In 
fact, the last of the above is not established nearly as 
firmly as the first three; nevertheless, we shall make use 
of it in developing our new formulation of the HFB+THO 
method. 

Beyond the point the HFB-I-HO solution fails to 
capture the physics of the coordinate-space results, es- 
pecially in the far asymptotic region. It is this incorrect 
large-r behavior that we now try to cure by introducing 
the THO basis. 



Approximation to the coordinate-space HFB 
local densities 



Our goal is to try to find an approximation to the ex- 
act (coordinate-space) HFB density that is based only on 
information contained in the HFB-f HO results. Towards 
that end, we make use of the WKB asymptotic solution 
of the single-particle Schrodinger equation for a given po- 
tential y(r)^ assuming that beyond the classical turning 
point only the state with the lowest decay constant k=2K 
contributes to the local density. Under this assumption, 
the logarithmic derivative of the density can be written 
as 



p'{r) 



p{r) 



V - 



1 V 



2k2 + V 



(3.2) 



where the first term comes from the three-dimensional 
volume element, while the next two correspond to the 
first- and second-order WKB solutions [s^l . The reduced 
potential V, 



n 



e{i+l) 2m Ze^ 



(3.3) 



is the sum of the nuclear, centrifugal, and Coulomb (for 
protons) contributions, with £ being the single-particle 
orbital angular momentum. 

In practical applications, it turns out that near Rm 
the next-to-lowest quasiparticle states still contribute to 
the local density p in a way that may be more important 
than the second-order WKB term shown in Eq. H3.2(l . 
Moreover, in deformed nuclei the quasiparticle states do 
not have good total angular momentum ^, so that several 
quasiparticles may contribute to the asymptotic density 
depending on their ^-content and the value of k. There- 
fore, we need a practical prescription to fix a reasonable 
approximate asymptotic form of the density with min- 
imal numerical effort but high reliability. This can be 
achieved by using in H3.2|l a reduced potential of the form 



V(r) 



C 2m Ze"^ 



(3.4) 



where the nuclear part V^v (which is small around and 
beyond Rm) is neglected, and the coefficient C is allowed 
to differ from its centrifugal barrier value £{i +1). The 
actual value of C is fixed by the requirement that the 
logarithmic derivative (|3.2|l coincides at the mid point 
Rm with the £—0 component of the HFB+HO density, 
i.e., with 



p{r) 



■it/2 



p(r, 6i)Pfco(cos(6i)) sin(6i)d6'. 



(3.5) 



Next, in order to make a smooth transition from the 
HFB+HO density p{r) in the inner region to the approx- 
imate asymptotic expression H3.2|) in the outer region, we 
introduce the following approximation p for the logarith- 
mic density derivative: 



P'jr) 
p{r) 



p(r) 

r. (^•." 



-r) 



for R„ 



for r < i?n 



-|-2V^ 



for r > i?max- 

(3.6) 

The coefficients a and 6, and the power s, are determined 
from the condition that the logarithmic derivative H3.6(l 
and its first derivative are smooth functions at the points 
i?min and i?max- Note that the first derivative of (|3.6() at 
-Rniin is automatically equal to zero, so that there is 
no need to introduce a fourth parameter to satisfy this 
condition. 

Having determined the smooth expression for the log- 
arithmic derivative of /5(r), we can derive the approxi- 
mate local density distribution p{r) by simply integrating 
Eq. lESJl. The result is 



p{r) 



for r < R 



A, 



-hr 



exp 



rRt 



B 



cxp[-2 f VK.'^+Vdr] 



2-5- ' 1-s 
for i?min <r < Rn 

for r > R„ 



(3.7) 

where the integration constants A and B are determined 
from the matching conditions for the density at points 
_Rniin and i?max, respectively. Finally, p(r) is normalized 
to the appropriate particle number. 

The approximate density 1)3. 7|) works fairly well for all 
nuclei that we have considered. This is illustrated in 
Fig. n where the approximate density p (circles) is seen to 
be in perfect agreement with the coordinate-space HFB 
results. 

It should be stressed that the above procedure is ap- 
plicable only when the number of shells is large enough 
that the HFB+HO density has a minimum at the point 
-Rmin- The minimum value of A^sh required to satisfy this 
condition depends on the particular deformations or on 
the nuclei considered. For the number of shells N^h ~ 20 
used in our calculations, the above condition is always 
satisfied. 
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C. LST function for HFB+THO calculations 

The starting point of our new and improved 
HFB+THO procedure is, thus, to carry out a standard 
HFB+HO calculation for the nucleus of interest, thereby 
generating its local density and its local £=0 density p{r) 
(|3.5|) . and then to use the method outlined in the previ- 
ous subsection to correct that density at large distances 
(see Eq. H3.7fl '). by calculating p{r). The next step is to 
define the LST 'lo] so that it transforms the HFB+HO 
£—0 density 1)3. 5|l into the corrected density of Eq. 1)3.7)1 . 
This requirement leads to the following first-order differ- 
ential equation. 



p{r) = 



p{n) dfjn) 
7^2 dn 



(3.8) 



which for the initial condition /(O) = can always be 
solved for f{TZ). 

Once the LST function has been so obtained, we need 
simply diagonalize the HFB matrices in the correspond- 
ing THO basis. Most importantly, no information is re- 
quired to build the THO basis beyond the results of a 
standard HFB+HO calculation. Since no further param- 
eters enter, there is no need to minimize the HFB+THO 
total energy. As a consequence, with this new methodol- 
ogy we are able to systematically treat large sets of nuclei 
within a single calculation. 

Despite the fact that the new HFB+THO method is 
simpler to implement than the earlier version, there are 
no discernible differences between the results obtained 
with the two distinct treatments of the LST function. 
Most importantly, the current formulation leads to the 
same excellent reproduction of coordinate-space results 
as did the previous one 0, 0| . 



IV. RESULTS 

In this section, we present the results of calculations 
performed for all particle-bound even-even nuclei with 
Z<108 and iV<188. The THO basis was implemented 
according to the prescription developed in the previous 
section. The k value used in the procedure was obtained 
in the following way. From the starting HFB+HO cal- 
culation, we determined k values separately for neutrons 
and protons, using Eq. ()3.1I) . We then associated the k 
value for the transformation with the smaller of kp and 
fc„. In this way, the THO basis is always adapted to the 
less-bound type of particle. The calculations were per- 
formed by building THO basis states from spherical HO 
bases with A'sh=20 HO shells and with oscillator frequen- 
cies of hujo = 1-2 X 41 MeV /A^/^. 

In order to meaningfully test predictions of nuclear 
masses for neutron-rich nuclei, we used the SLy4 Skyrme 
force parameterization [23j, as this was adjusted with 
special emphasis on the properties of neutron matter. At 
present, there also exist Skyrme forces that were adjusted 
exclusively to nuclear masses |35| . These forces were 



used within a calculation scheme that was not focused on 
weakly-bound nuclei. In the pairing channel, we used a 
pure volume contact pairing force V^{r, r') = Vb(5(r — r') 
with strength Vo=— 167.35 MeV fm'^ and acting within 
a phase space limited by a cut-off parameter 19j of 
emax=60MeV. 

Figure |5] summarizes the systematic results of our cal- 
culations, both for ground state quadrupole deformations 
(upper panel) and for two-neutron separation energies 
(lower panel). For this figure, calculations for a given 
mass number A were carried out for increasing (decreas- 
ing) N—Z, up to the nucleus with positive neutron (pro- 
ton) Fermi energy. Furthermore, for each nuclide, three 
independent sets of HFB+THO+LN calculations were 
performed, for initial wave functions corresponding to 
oblate, spherical, and prolate shapes, respectively. De- 
pending on properties of a given nucleus, we could there- 
fore obtain one, two, or three solutions with different 
shapes. For each obtained solution we performed a PNP 
calculation of the total energy. The lowest of these en- 
ergies for a given nucleus was then identified with the 
ground-state solution. 

Calculations of a microscopic mass table are greatly 
helped by taking advantage of parallel computing. We 
have used two IBM-SP computers at ORNL: Eagle, a 
1 Tflop machine, and Cheetah, a 4 Tflop machine (1 
Tflop = 1 X 10^^ operations/second). The code per- 
forms at 350 Mflop/processor on Eagle. We created a 
simple load-balancing routine that allows us to scale the 
problem to 200 processors. We are able to calculate the 
entire deformed even-even mass table in a single 24 wall- 
clock hour run (or approximately 4,800 processor hours). 
A complete calculated mass table is available online in 
Ref. [ig. 

The ground-state quadrupole deformations /? dis- 
played in Fig. 12 (upper panel) were estimated from 
the HFB+THO+LN total quadrupole moments and rms 
radii through a simple first-order expression 's^l . In that 
panel, all even-even nuclei with negative Fermi energies, 
A„<0 and Ap<0, are shown. In the lower panel, showing 
two- neutron separation energies S2n, results are shown 
for those N and Z values for which the nuclides with 
both N and N —2 have A„<0. Note that on the proton- 
rich side the lighter of them may have Ap>0; nevertheless, 
we show these points to make the proton drip line in the 
S2n panel identical to that of the quadrupole deforma- 
tion panel. Of course, on the proton drip line values of 
S2n are large and not very illuminating. 

Table summarizes our results for even-even nuclei 
along the two-particle drip lines. More specifically, for 
each value of Z, the results for the lightest isotope with 
Ap<0, and the heaviest isotope with A„<0 are presented. 

As can be seen from Fig.Eland TablelH our calculations 
produce several particle-bound even-even nuclei (i.e., nu- 
clei with negative Fermi energies) that at the same time 
have negative two-proton (or two-neutron) separation en- 
ergies. Such an effect was already noticed in light nuclei 
m Ref. M- 

The current calculations suggest it may be 
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FIG. 2: Quadrupole deformations /3 (upper panel) and two-neutron separation energies S2n in MeV (lower panel) of particle- 
bound even-even nuclei calculated within the HFB-f THO method with Lipkin-Nogami correction followed by exact particle 
number projection. The Skyrme SLy4 interaction and volume contact pairing were used. 



generic, occurring near both the two-neutron and two- 
proton drip lines and for nuclei as light as ^^Mg and 
as heavy as ^^^Dy. It seems to be related to the fact 
that the Fermi energies pertain to stability with respect 
to particle emission of a given configuration or shape, 
namely that of the ground state. In many of the cases 
in which we observe this phenomenon, (a) the neighbor- 
ing even-even nucleus, the one to which it would decay 
by two-nucleon emission, has two distinct shapes, each 



with negative Fermi energies, (b) the ground state of that 
neighboring nucleus has a shape that is different than 
that of the parent nucleus, (c) the shape of the excited 
bound configuration is the same as that of the parent 
nucleus, and (d) decay to the excited configuration is en- 
ergetically forbidden. 

The precise results of course depend sensitively on 
properties of the interaction, both in the particle-hole 
and particle-particle channels. Despite its many good 
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TABLE I: Results of the HFB+THO calculations for drip-line nuclei with the SLy4 Skyrme force and volume delta pairing 
force. The left and right columns show results for proton and neutron drip-line isotopes from He to Pb. For both drip lines we 
show deformations (3, Fermi energies A (in MeV), two-particle separation energies (in MeV), and neutron and proton pairing 
gaps (in MeV). 
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features, the force we use is far from perfect. For exam- 
ple, the positions we obtain for the two-neutron drip Unes 
in the Be and O isotopes are not correct. In addition, the 
method itself has limitations, as it leaves out potentially 
important effects beyond mean field. Despite these limi- 
tations, we feel it is nevertheless worthwhile to point out 
some of the interesting new physical situations that are 
predicted in these calculations and which may therefore 
occur in weakly-bound systems. The above example of 
nuclei that are formally beyond one of the two-particle 



drip lines but nevertheless are localized and do not spon- 
taneously spill off a nucleon is just one of several. We 
will now discuss in greater detail some specific isotopic 
chains to see how this and other interesting exotic new 
features emerge. 

We focus our discussion on the heaviest isotopes of 
four isotopic chains; neon, magnesium, sulfur, and zinc 
(see Figs.lHHOl respectively). The figures show the Fermi 
energies A^^, A^, and [see Eqs. ()2.13f) and (|2.14|) ]. 
and the total binding energies, obtained in constrained 
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FIG. 3: Neutron Fermi energies A (upper panels) and the total 
binding energies (lower panels) calculated for ^"Ne, ^^Ne, and 
^''Ne as functions of the quadrupole deformation f3. 



HFB+THO+LN+PNP calculations as functions of the 
quadrupole deformation /3 for the last three particle- 
bound isotopes of the respective chains. In each figure, 
the binding energies of the last three isotopes are shown 
on a common energy scale. As a reminder, two neu- 
tron separation energies can be readily obtained from 
the binding energies according to S2n ~ E{Z, N ~ 2) — 
E{Z,N). 

We should note that the minima of the constrained 
energies need not exactly correspond to the PNP of the 
HFB-t-THO-)-LN minima, which were used in Fig. Inland 
Tabled Indeed, in the constrained calculations the de- 
formation serves as an additional variational parameter 
for the variation after PNP. Optimally, the full variation 
after projection should be performed, which, however, 
requires a much larger numerical effort, and is left to fu- 
ture work. Such an optimal method will also remove the 
ambiguities related to the definition of the Fermi energy, 
discussed in Sec.m At present, we illustrate these ambi- 
guities by showing in Figs. I^HHI the three possible values 
of the Fermi energy, A„, A^, and A". 

Consider first the Ne isotopes, for which the results are 
shown in Fig. O For the SLy4 interaction that we use, a 
strong shell gap at iV=20 persists up to the heaviest iso- 
topes of Ne, and this produces a stiff spherical minimum 
for '^"Ne. Adding two neutrons gives rise to the nucleus 
■^^Ne, which is particle-bound (A„<0), but at the same 
time two-neutron unstable (<S'2n<0). [Note that this nu- 
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FIG. 4: Same as in Fig.^lbut for ^**Mg, ''"Mg, and ''^Mg. 



cleus does not exactly fit into the picture given earlier for 
such nuclei.] Interestingly, when we add two more neu- 
trons, we obtain a strongly (prolate) deformed particle- 
bound ground configuration in '^^Ne, which is again two- 
neutron stable (S'2n>0). 

Next we turn to the Mg isotopes, for which results 
are presented in Fig. 0] In ""^Mg the neutron Fermi en- 
ergies A„ have negative values for all deformations, so 
that the configurations for all deformations are particle- 
bound, with the prolate minimum being slightly lowest. 
The same is also true for the next nucleus ''^Mg where the 
ground state deformation changes from prolate to oblate. 
It is clear from Fig.^jthat in ''^Mg the two-neutron sepa- 
ration energy is negative; however, since ^^Mg and '^''Mg 
have different shapes in their ground states, the real pro- 
cess of emitting two neutrons may occur towards the 
shape isomer in ^"^Mg. (The situation will be even more 
complicated if the oblate minimum in ''^Mg is unstable 
to triaxial deformations, i.e., it is a saddle point.) 

The results for the S isotopes are given in Fig. |5l 
Here, the spherical HFB-t-THO+LN minimum in '"^^S is 
shifted in the constrained PNP calculations towards a 
small oblate deformation. All shapes appear to be very 
weakly particle-bound, and have negative two-neutron 
separation energies at the same time. It is obvious that 
in the case of so poorly defined a minimum, its pre- 
cise location is not relevant and full configuration mix- 
in g , e.g., w ithin the generator coordinate method (GCM) 
|30ll37ll3q . should be applied. This complication specific 
to weakly-bound nuclei is related to the fact that it is not 
clear how to take into account in the GCM the regions 
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FIG. 5: Same as in Fig.Olbut for ^^S, ^°S, and ^^S. 



of the collective coordinate corresponding to A>0, hence 
to particle-unbound states. 

In the results for the zinc isotopes (Fig. (BJ, we see 
strong competition between oblate, prolate, and spheri- 
cal shapes. In ^^Zn, all shapes are particle-bound and 
the ground state is oblate. The situation changes in 
^^Zn, where the oblate configuration, though lowest in 
energy, becomes particle-unbound and the prolate min- 
imum becomes the ground-state configuration. Though 
this ground state is two-neutron unstable (5*271 <0), its 
decay to the ground state of ^^Zn may be hindered by 
the shape change. Finally, in ^'^"Zn the particle-stable 
prolate ground state is also two-neutron unstable. Hence 
in this isotopic chain the last two even isotopes are un- 
stable with respect to two-neutron emission. 

In heavier nuclei near the neutron drip line, we often 
obtain particle-stable and two-neutron-unstable isotopes 
right after closed neutron magic shells. As in Ne, this 
reflects the fact that strong shell gaps persist up to the 
heaviest isotopes in a chain when the calculations are 
based on the SLy4 interaction. In the N—126 isotopes 
of Ce and Nd, for example, the ground-state configura- 
tions are strongly spherical. In the neighboring iV=128 
isotopes, these spherical configurations become particle 
unbound. However, in these same isotopes, there are 
strongly prolate particle-bound configurations with very 
large negative two-neutron separation energies (see Ta- 
bleQ}. An analogous situation occurs in the iV=186 and 
188 drip-line nuclei, where the last two even isotopes may 
have particle-bound prolate states with unbound spheri- 
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cal configurations. 

Strong SLy4 neutron magic numbers also result in the 
characteristic non-monotonic behavior of the S2n values 
(Fig. 12). Indeed, lines of constant S2n often follow de- 
creasing Z with increasing N, which is particularly con- 
spicuous near 7V=126. This effect even creates a small 
peninsula of stability near iV=140. Such strong neutron 
closed shells could create the well-known deficiencies in 
the r-process abundances [Sflj ]. 



V. CONCLUDING REMARKS 

In this paper, we have reported the development of an 
improved version of the configuration-space HFB method 
expanded in a transformed harmonic oscillator basis. In 
its current form, the method can be used reliably in sys- 
tematic studies of wide ranges of nuclei, both spherical 
and axially deformed, extending all the way out to the 
nucleon drip lines. The key step was the development 
of a prescription for choosing a reliable transformation 
function to define the THO basis that does not require 
variational optimization. The current prescription only 
involves information from a preliminary configuration- 
space HFB calculation carried out in a harmonic oscilla- 
tor basis. The transformation function is then tailored 
to correct the asymptotic properties of the HFB-I-HO re- 
sults. The resulting HFB-I-THO theory accurately re- 
produces results of coordinate-space HFB theory, where 
available, and also reproduces the results obtained with 
an earlier version of the transformation that had to be 
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optimized separately for each nucleus. 

As a first application of the new HFB+THO method- 
ology, we carried out a systematic study of all even- 
even nuclei having Z<108 and A^<188. Variation after 
particle-number projection was approximately included 
using the Lipkin-Nogami method, with exact projection 
performed for the final self-consistent solutions. We fo- 
cussed our discussion on those nuclei that are very near 
the nucleon drip lines, finding that in several regions 
of the periodic table there exist nuclei that are stable 
against one-particle emission but unstable against pair 
emission. We showed that invariably this is associated 
with a shape change in the ground state. For example, 
while two-particle emission to the configuration of the 
daughter with the same shape as the parent is forbid- 
den, a decay to the ground state having a different shape 
can nevertheless occur. The associated change in shape 
may conceivably lead to sufficient hindrance of the decay, 
hence the longer lifetime. Consequently, it is conceivable 
that there exist nuclei that formally live beyond the neu- 
tron drip line but can be observed experimentally. This 
phenomenon, which had earlier been noted in calcula- 
tions of light nuclei, is now seen to be a more common 
feature of nuclei near the neutron drip line. 

In the description of very weakly-bound systems, small 
changes in the results can have important consequences, 
changing for example the precise locations of the drip 
lines. It is important, therefore, to continue to improve 
the current HFB+THO methodology to accommodate 



effects not presently being included. Particularly impor- 
tant could be effects that arise beyond mean field. It is 
also important to develop the new HFB-I-THO formalism 
for application to odd-mass systems, including the effects 
of Pauli blocking. But most crucial, in our opinion, is to 
develop new-generation energy density functionals that 
will allow for more reliable predictions of the properties 
of exotic nuclei. Work along these various lines is cur- 
rently underway. 
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